ConvertTransverseMercatorToGeodetic Subroutine

SUBROUTINE ConvertTransverseMercatorToGeodetic &
!
(x, y, k, centM, lat0, a, e, eb, falseN, falseE, lon, lat)

USE Units, ONLY: &
!Imported parametes:
pi

USE StringManipulation, ONLY: &
!Imported routines:
ToString

IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT (IN) :: x !!easting coordinate [m]
REAL (KIND = float), INTENT (IN) :: y !!northing coordinate [m]
REAL (KIND = float), INTENT (IN) :: k !!scale factor
REAL (KIND = float), INTENT (IN) :: centM !!central meridian [radians]
REAL (KIND = float), INTENT (IN) :: lat0 !!latitude of origin [radians]
REAL (KIND = float), INTENT (IN) :: a !! semimajor axis [m]
REAL (KIND = float), INTENT (IN) :: e !! eccentricity
REAL (KIND = float), INTENT (IN) :: eb !! second eccentricity
REAL (KIND = float), INTENT (IN) :: falseN !!false northing
REAL (KIND = float), INTENT (IN) :: falseE !!false easting


!Arguments with intent (out):
REAL (KIND = float), INTENT (OUT) :: lon !!geodetic longitude [radians]
REAL (KIND = float), INTENT (OUT) :: lat !!geodetic latitude [radians]

!Local variables:
REAL (KIND = double) ::  c  !!Cosine of latitude 
REAL (KIND = double) :: es !!eccentricity squared
REAL (KIND = double) :: ebs !!second eccentricity squared
REAL (KIND = double) :: de  !!Delta easting - Difference in Easting (Easting-Fe)
REAL (KIND = double) :: dlam !!Delta longitude - Difference in Longitude
REAL (KIND = double) :: eta !!constant - ebs *c *c 
REAL (KIND = double) :: eta2 
REAL (KIND = double) :: eta3
REAL (KIND = double) :: eta4
REAL (KIND = float) :: ftphi !! Footpoint latitude
INTEGER (KIND = short) :: i  !! Loop iterator
REAL (KIND = double) :: s !!Sine of latitude
REAL (KIND = double) :: sn !!Radius of curvature in the prime vertical
REAL (KIND = double) :: sr !!Radius of curvature in the meridian
REAL (KIND = double) :: t !!Tangent of latitude
REAL (KIND = double) :: tan2
REAL (KIND = double) :: tan4
REAL (KIND = double) :: t10 !!Term in coordinate conversion formula
REAL (KIND = double) :: t11 !!Term in coordinate conversion formula
REAL (KIND = double) :: t12 !!Term in coordinate conversion formula
REAL (KIND = double) :: t13 !!Term in coordinate conversion formula
REAL (KIND = double) :: t14 !!Term in coordinate conversion formula
REAL (KIND = double) :: t15 !!Term in coordinate conversion formula
REAL (KIND = double) :: t16 !!Term in coordinate conversion formula
REAL (KIND = double) :: t17 !!Term in coordinate conversion formula
REAL (KIND = double) :: tmd !!True Meridional distance 
REAL (KIND = double) :: tmdo !!True Meridional distance for latitude of origin 

! Local parameters:
!! Maximum variance for easting and northing values for WGS 84.
REAL (KIND = float), PARAMETER :: deltaEasting = 40000000.0
REAL (KIND = float), PARAMETER :: deltaNorthing = 40000000.0

!------------end of declaration------------------------------------------------


! Check consistency of input coordinates
IF ( x < (falseE - deltaEasting) .OR. &
     x > (falseE + deltaEasting)      ) THEN
 
   CALL Catch ('error', 'GeoLib',   &
			   'Converting Transverse Mercator to Geodetic: &
			    easting out of range' ,  &
			    code = consistencyError, argument = ToString(x) )
END IF

IF ( y < (falseN - deltaNorthing) .OR. &
     y > (falseN + deltaNorthing)      ) THEN
 
   CALL Catch ('error', 'GeoLib',   &
			   'Converting Transverse Mercator to Geodetic: &
			    northing out of range' ,  &
			    code = consistencyError, argument = ToString(y) )
END IF

! True Meridional Distance for latitude of origin 
tmdo = MeridionalDistance (lat0, a, e)

!Origin
tmd = tmdo +  (y - falseN) / k 

!First Estimate
sr = Rho (0., a, e)
ftphi = tmd / sr

DO i = 1,5
  t10 = MeridionalDistance (ftphi, a, e)
  sr = Rho (ftphi, a, e)
  ftphi = ftphi + (tmd - t10) / sr
END DO

! Radius of Curvature in the meridian
sr = Rho (ftphi, a, e)

! Radius of curvature in the prime vertical
sn = Nu (ftphi, a, e)

!Eccentricities squared
es = e**2.
ebs = eb**2.

! Sine Cosine terms
s = SIN(ftphi)
c = COS(ftphi)

! Tangent Value
t = TAN(ftphi)
tan2 = t * t
tan4 = tan2 * tan2
eta = ebs * c**2.
eta2 = eta * eta
eta3 = eta2 * eta
eta4 = eta3 * eta
de = x - falseE
IF ( ABS(de) < 0.0001) THEN
   de = 0.0
END IF

! Latitude 
t10 = t / (2. * sr * sn * k**2.)
t11 = t * (5.  + 3. * tan2 + eta - 4. * eta**2. - 9. * tan2 * eta) / &
          (24. * sr * sn**3. * k**4.)
t12 = t * (61. + 90. * tan2 + 46. * eta + 45. * tan4 &
           - 252. * tan2 * eta  - 3. * eta2 + 100.   & 
           * eta3 - 66. * tan2 * eta2 - 90. * tan4   &
           * eta + 88. * eta4 + 225. * tan4 * eta2   &
           + 84. * tan2* eta3 - 192. * tan2 * eta4 ) &
           / ( 720. * sr * sn**5. * k**6. )
t13 = t * (1385. + 3633. * tan2 + 4095. * tan4 + 1575. * t**6. ) & 
           / (40320. * sr * sn**7. * k**8. )

lat = ftphi - de**2. * t10 + de**4. * t11 - de**6. * t12 + de**8. * t13

t14 = 1. / (sn * c * k)

t15 = (1. + 2. * tan2 + eta) / (6. * sn**3. * c * k**3.)

t16 = (5. + 6. * eta + 28. * tan2 - 3. * eta2 &
       + 8. * tan2 * eta + 24. * tan4 - 4.    & 
       * eta3 + 4. * tan2 * eta2 + 24.        & 
       * tan2 * eta3) / (120. * sn**5. * c * k**5. )

t17 = (61. +  662. * tan2 + 1320. * tan4 + 720. & 
       * t**6. ) / (5040. * sn**7. * c *k**7. )
       
!Difference in Longitude
dlam = de * t14 - de**3. * t15 + de**5. * t16 - de**7. * t17

!Longitude
lon = centM + dlam

! Check errors
IF ( ABS (lat) > 90. * degToRad ) THEN
   CALL Catch ('error', 'GeoLib',   &
    'Converting Transverse Mercator to Geodetic: &
    latitude out of range' ,  &
    code = consistencyError, argument = ToString(lat*radToDeg)//' deg' )
END IF
   
IF( lon > pi ) THEN
  lon = lon - (2. * pi)  
       
ELSE IF (lon < -pi ) THEN
  lon = lon + (2. * pi)
END IF

IF( ABS (lon) > pi ) THEN
     CALL Catch ('error', 'GeoLib',   &
     'Converting Transverse Mercator to Geodetic: &
     longitude out of range' ,  &
     code = consistencyError, argument = ToString(lon*radToDeg)//' deg' )
END IF

! Check for distortion and send warning.
!Distortion will result if Longitude is more than 9 degrees 
!from the Central Meridian at the equator
!and decreases to 0 degrees at the poles
!As you move towards the poles, distortion will become more significant 
IF ( ABS(dlam) > (9.0 * degToRad) * COS(lat) ) THEN
   CALL Catch ('warning', 'GeoLib',   &
     'Converting Transverse Mercator to Geodetic: &
     possible distortion in longitude computation ' ,  &
      argument = ToString(lon*radToDeg)//' deg' ) 
END IF

END SUBROUTINE ConvertTransverseMercatorToGeodetic